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Abstract 

We present a qualitative diagnostic based on the continuous wavelet transform, for 
the detection of irregular behaviour in time series of particle simulations. We apply 
the method to three qualitatively different gravitational 3-body encounters. The 
intrinsic irregular behaviour of these encounters is well reproduced by the presented 
method, and we show that the method accurately identifies the irregular regime in 
these encounters. We also provide an instantaneous quantification for the degree of 
irregularity in these simulations. Furthermore we demonstrate how the method can 
be used to analyze larger systems by applying it to simulations with 100-particles. 
It turns out that the number of stars on irregular orbits is systematically larger for 
clusters in which all stars have the same mass compared to a multi-mass system. The 
proposed method provides a quick and sufficiently accurate diagnostic for identifying 
stars on irregular orbits in large scale A^-body simulations. 
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1 Introduction 



Gravitational A^-body systems are intrinsically chaotic, at least for N > 2 [1]. 
For N = 2 Newtonian systems are regular, but post Newtonian dynamics [2, 3, 4] 
can already reveal chaotic behaviour [5,6] by = 2. Chaoticity in a self grav- 
itating systems with N > 2 can reveal itself on a very short time scales, of the 
order of orbital periods. On the other hand, in some systems like the plane- 
tary system around the Sun, chaoticity only reveals itself on a time scale of 
billions of orbits [7,8]. In a stellar cluster which does not contain a dominant 
central mass, orbits can be chaotic on a much smaller time scale. In such an 
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environment seemingly regular orbits can suddenly become highly irregular, 
to return later to regular again [9,10]. 



In this study we focus on the characterisation of orbits which show irregular 
behaviour on a short (dynamical) time scale, and not on a long (relaxation) 
time scale. Star clusters, with contain between about 100 and 0(10''') stars, 
are governed by microscopic few body interactions [11]. Analysis of such sys- 
tems is often hampered by this internal microscopic physics and a qualitative 
indicator of chaotic behaviour could enormously assist in the understanding of 
large scale gravitational A'"-body simulations. In particular, since it is thought 
that stars on irregular orbits in self gravitating iV-body systems have an im- 
portant effect of the bulk properties of such systems [12]. 



We present a transparent diagnostic to qualify chaotic behaviour of individual 
trajectories in A^-body systems. In addition, the method has some quantita- 
tive qualities. The application of this method ranges from 3-body interactions 
to simulations of entire star clusters and galaxies. For clarity, we define ir- 
regular orbits as orbits which are deterministic though sensitive to the initial 
conditions and which cannot be described as a sum of periodic motions. 



The motion of a star in an A-body simulation can change from regular to 
irregular (and vice versa) due to gravitational interactions with other stars. 
Irregularity of an orbit then is a local quantity, and this forces us to devi- 
ate from the conventional methods based on Lyapunov [13] numbers such as 
explained in [14,15,16]. In addition, the calculation of Lyapunov exponents 
is costly, may require reruns and are therefore less suited for a direct diag- 
nostic, whereas we are predominantly interested in a diagnostic that can be 
evaluated at runtime. The Lyapunov indicator, however, is well suited for dis- 
criminating between ordered and weak chaotic motion, like planetary systems, 
whereas we are in particular interested in very chaotic systems [14]. Others 
methods for detecting chaos in Hamiltonian systems, such as SALI [17,18], 
Fourier Transform [19,20,21], Poincare Section [22], the zero-one method [23] 
and the geometric indicator [24] are also less suited for analysing gravitational 
A^-body simulations during run time, since they also are not practical in pro- 
viding an instantaneous quantification of the chaoticity of the system we are 
interested in here. 



In §2 we discuss the basics of the Continuous Wavelet Transform (CWT), 
which we present in a package called CWaT. We then proceed by applying 
this method to several A-body simulations in § 3, to conclude in § 4. 
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2 Method for detecting chaos 



2.1 The continuous wavelet transform 



The continuous wavelet transform [25] (CWT) gives the time-frequency rep- 
resentation of a time series /(t) by fitting a wavelet ^ to it at subsequent 
points in time t. We have used the Morlet-Grossmann wavelet [26], which is 
described by 

= e^'"°*e-*'/''^'. (1) 



Here cr is a measure of the spread in time and loq the base frequency of the 
wavelet. The wavelet is a periodic sinusoidal signal enveloped by a Gaussian, 
from which we can construct a family of wavelets using a scaling factor a and 
a time shift /3 

- f^") , Pen, a>0. (2) 



We define the CWT as 

oo 

T(a,(3)^ I f{t)^^^,p){t)dt. (3) 



Here is the complex conjugate of ^ and f{t) a time scries. Fitting the 
wavelet to f{t) at discrete moments in time {(3 G [0,t]) yields T{a,t). Wc 
express T{a,t) in terms of frequency uj rather than scale a (since we are more 
accustomed to the former) using the relation 

uj^—. (4) 
a 

This leads to a time-frequency representation T{u, t). For an extensive discus- 
sion see [27] . 



2.2 Analysing the CWT 



The main features of the time series can be described by its instantaneous 
frequencies u}{t). These can be extracted from T{u},t) by selecting its local 
maxima. We only consider all local maxima larger than a threshold r = 0.95 
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of the global maximum. In this we take a different approach than [28,29] who 
explicitly consider the global maximum. Using a lower threshold enables us to 
filter most of the noise from the CWT without losing the relevant information 
on which the x can be determined. We also tried r = 0.93 and r = 0.97, which 
gave very similar results in the resulting value of x- 

We find the extrema by simulated annealing [30,31], following [32]. Connecting 
the local maxima results in a number of curves in the time- frequency domain, 
which we call ridges. The time series is subsequently analysed by studying the 
ridges, enabling a qualitative analysis of the behaviour of the time series. This 
is considerably more practical than using the CWT directly, an example is 
presented in fig. 2. 

A periodic time series (i.e. with a constant frequency) results in a uninter- 
rupted horizontal ridge indicating the fundamental frequency P and possibly 
overtones spaced at periods P/n where n is integer. A quasi-periodic time- 
series can be described by a sum of periodic time series, each with its own 
ridge at a specific frequency. Such a time series is then represented by multi- 
ple horizontal ridges. Irregular time series on the other hand lead to curved 
ridges of limited length [33] . 



2.3 A quantitative measure of chaos 

To bring the chaos detecting qualities of the method one step further we 
quantify the number of ridges and their curvature in the irregularity indicator 
X- We define x at any moment in time by: 

I n{t) 

xW = ;^gW. (5) 



Here n(t) is the number of ridges at time t and Sj a measure for the slope of 
ridge i. Using a least squares fit we calculate the slope of the ridge. Because 
X would be dominated by large values of the slope we use an empirical cutoff. 
The slope is normalised to unity to assure that Sj G [0, 1]. 

The ridges of periodic and quasi-periodic time series are represented by hor- 
izontal ridges, they imply ~ (see §2.2). We note however that the 
proposed indicator (see Eq. 5), cannot give a yes/no definition about the oc- 
currence of local exponential instability, but it is capable of giving a measure 
of the degree of instability. 
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3 Application 



To demonstrate the effectiveness of CWaT we apply it to three gravitational 
3-body interactions. Such interactions are irregular [34], with some regular 
characteristics (see [35,36] for an interesting earlier analysis of this sort). 

In each of the three following examples we inject a single object (star) in a sys- 
tem consisting of two objects (stars) which are bound in a binary. A time series 
is constructed by projecting the orbit of each of the stars on three orthogonal 
axes (see also [19]). The stellar orbits are computed using the starlab soft- 
ware environment [37,38] which is publicly available from 
http://www.ids.ias.edu/~starlab This A-body integration package cal- 
culates the orbits of the stars with individual time steps [39] using a fourth- 
order Hermite predictor-corrector scheme [40]. A regular time series is con- 
structed by Hermite interpolation at evenly spaced time intervals. The time 
step adopted for CWT is 1/1024*"^ of the total integration time. Further anal- 
ysis is carried out with the S-WAVE software package [32], to ultimately cal- 
culate the CWT and the ridges of the time series. 

In the following two §§ we discuss the chaotic behaviour in one non-resonant 
preserving encounter and one resonant encounter. In the first encounter the 
incoming star escapes after a close encounter with a binary. In the following 
example, in §3.2 (see [41] for more examples, and even some very peculiar or- 
bits in [42]), the incoming star has multiple close interactions with the binary, 
we call such interactions 'resonant'. The example shown is a democratic reso- 
nance which is characterised by all three stars interacting strongly with mul- 
tiple close passages; each star generally interacts with comparable frequency 
with any of the other stars. 

We use the following convention: star 7^1 and 7^2 form the initial binary 
which is encountered by star 7^3. The time series give the position of the star 
with respect to the centre of mass of the 3-body system. The units are given 
in standard A-body units [43]. 

3. 1 Irregular motion in a non-resonant gravitational interaction 

In the first example we select an interaction between a binary, consisting of 
star 7^1 and jj^2 to interact with a single incoming star, 7^3. After inter- 
acting with the binary, star 7^3 escapes without disrupting the binary. The 
incomer, however, may change the binding energy of the binary. Precise initial 
conditions are presented in Tab. .1. 

In fig. . 1 we present a schematic diagram of the interaction. In this example the 
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binary has a close encounter with the incomer at t ~ 150. After the encounter 
the incoming star is deflected along a hyperbolic path. 

The analysis using CWaT of one of the nine (three stars and three axes) time 
series is given in fig. 2. The top panel shows the time series of the position 
of star projected on the y-axis in A''-body spacial coordinates versus A^- 
body time units. In the first 150 A^-body time units, t = [0, 150], the binary 
moves upward (in y) with respect to the centre of mass because the incomer 
(star #3) is approaching it from above. The projection of the circular motion 
results in the sinusoidal behaviour of the time series. The motion of the binary 
with respect to the centre of mass is downwards after the encounter with the 
incomer. 

The second panel of fig. 2 shows T{u;,t). Its most prominent feature is two 
arch-like shapes in the period region 10 to 100 spanning t = [30, 150] and 
t = [150, 330]. These arches are the sum of one arch spanning t = [30, 330] due 
to boundary effects [44] and a triangular feature centred at t ^ 150 due to the 
centre of mass motion of the binary. These regions of high intensity result in 
ridges with periods P ~ 86 and P ~ 61 in the third panel. 

The ridge at P ~ 6.3 in the third panel indicates the period of the binary 
system, which after t ~ 150 has reduced to P ~ 6.0. This ridge is joined by 
another on the interval t — [150, 330] at a shorter period of P ~ 3. The lower 
period ridge is an overtone caused by an increase of the eccentricity of the 
binary from e = to e ~ 0.38. 

As we discussed in Section 2.3, the bottom panel of fig. 2, shows that % ~ 
during the entire interaction in accordance with the periodic behaviour of the 
system. The short fiy-by interaction has not let to a qualitative indication of 
irregular behaviour in the system. 

A similar analysis can be performed for any star and for each projection. 
The result will be quite similar, supporting the argument that the presented 
method is independent of the selection of a fundamental plane. 

3.2 Irregular motion in a democratic resonant interaction 

In a democratic resonance encounter the incoming star interacts closely with 
both components of the initial binary, until one of the stars eventually escapes. 
The initial conditions used for this example are presented in Tab. 2 and a 
schematic representation of the interaction is shown in fig. 3, and the analyse 
using CWaT is presented in fig. 4. 

This interaction is considerably more complicated that the previous two ex- 
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Fig. 1. Schematic representation of a non-resonant three body interaction. The 
interacting stars on the left are indicated by 1, 2 and 3, each of which progresses in 
time along the solid curve to the right. Initially star 1 and 2 are close together, as 
is also indicated in the representation. A close encounter with star 3 sometime half 
way the encounter is marked by the three lines coming close together. The label 
"123" identifies the moment of closest approach between the three stars. To the 
top of the schematic representation we plotted a bar with the letters 'A' through 
'C, which correspond to various regions of different behaviour. A more quantitative 
analysis is presented in fig. 2 



Fig. 2. Analysis of the non-resonant 3-body interaction. The schematic interaction 
is presented in fig. 1 and the initial conditions are provided in tab. 1. The top panel 
gives the time series y{t) of star #2. The sinusoidal line to the bottom right shows 
a projection of the orbital motion of star ^2. The second panel contains the CWaT 
for this time scries, followed by the ridge extraction in the third panel from above. 
The bottom panel shows the x according to equation 5. 

amples. In fig. 3 and 4 we see how the incomer, star ^3, is captured by the 
binary (consisting of star ^1 and ^2) at t ~ 400 and interacts several times 
strongly with the other two stars before being ejected again at t ~ 1700. The 
initial binary has a period P ~ 6.3, but after t ~ 400 the CWT becomes very 
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Fig. 3. Schematic representation for the democratic resonant 3-body interaction. 
The analysis of this interaction is given in fig. 4. 



Fig. 4. Analysis of the democratic resonant 3-body interaction from fig. 3. The 
analysis is based on the time-series x{t) of star #1. 

irregular. The most irregular regime is noticeable in the analysis for x in the 
bottom panel of fig. 4 at t = [450,1700], with a peak near t = [750,1000], 
where values of x ^ 0.6 are reached. 

The strong democratic resonance interaction in this § provides an excellent 
example to compare the CWaT method with Poincare sections, which we will 
do in the subsequent §. 

3.3 Poincare section of a democratic resonant interaction 

Poincare sections have, like the CWaT diagnostic, the advantage that a qual- 
itative impression of the irregularity in time series can be obtained instanta- 
neously, id est; at run time. The main disadvantage of Poincare sections is the 
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Fig. 5. Poincare section of the velocity components of the democratic resonant 
3-body interaction shown in fig. 3 and further analysed in fig. 4. Here we show 
Vx and Vz for star #1 at each moment that Vy becomes negative. The letters 'A', 
'B', 'C and 'E' correspond to the regions indicated in fig. 3. Scattered asterisks 
between region 'B' and 'C belong to time interval indicated with 'D'. The little 
squares indicate the short intermittent period between 'A' and 'B'. 



required choice of a fundamental plane. This hinders analysis of complicated 
orbits in larger iV-body systems, as there is no simple automated procedure 
to find a fundamental plane through which the Poincare section can be taken. 
This is in particular a problem in systems with hierarchical substructure, such 
as the motion of moons around planets which again move around a star. The 
CWaT method does not suffer from this complication and, in addition, pro- 
vides an unbiased indicator for the degree of irregularity. We here compare the 
CWaT method with Poincare sections to demonstrate that global behaviour 
is reproduced in both methods. 

Poincare sections arc compiitcd using the velocity components of a star along 
a fundamental plane and record the moments a star passes through that plane. 
A stable binary in this representation gives a linear relation between the veloc- 
ity components at each intersection. For the democratic resonant interaction 
of § 3.2 we use the velocity components of star #1 in the VxVz-p\.Sine at the mo- 
ments Vy becomes negative. The centre of mass of the 3-body system remains 
in the origin throughout the simulation. The results of this are presented in 
fig 5. 

In fig 5 we indicate the five epochs by the letters 'A' through 'E'. The same 
designation is used in fig. 3. The points is regions 'A' and 'E' are clustered 
together, indicating that the orbit of the star before ('A') and after ('E') is 
well behaved. During the resonant encounter ('B', 'C and 'D') each of these 
resonant epochs can be grouped with hindsight from the CWaT method. In 
fig. 5 we use different symbols to identify the various resonant epochs. 
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Fig. 6. Evolution of the radii which contains 5% (bottom thick sohd curves), 50% 
(middle solid curve) and 75% (thin solid curve) of the mass for simulation model 
EqM (all masses are equal). For model S, with a Salpctcr initial mass distribution 
for the individual particles, the same data is presented as the thick dotted curve 
(5% Lagrangian radius) and the two thin dotted curves representing the 50% and 
the 75% Lagrangian radii. 

3.4 Irregular motion in a 100-body system 

After having confirmed that the CWaT method enables a reliable detection 
of irregular orbits in small N-body systems it is time to focus on the larger 
systems. Here we will apply CWaT to self gravitating A^-body systems with 
100 particles. In particular we report here the result for several of such simula- 
tions. The initial realizations of the particle clusters ware generated as follows. 
We make a distinction between two basic models, model EqM and model S. 
Model EqM was constructed by distributing 100 objects of the same mass in a 
Plummer [45] Spherical density distribution. The velocity of all point masses 
are subsequently scaled such that the entire system is in virial equilibrium, 
and to A'-body units [43]. We run the system for 100 A^-body time units, 
which corresponds to about 35 dynamical time units. 

In fig. 6 we present the evolution of several fundamental radii for simulation 
model EqM which contain 5%, 50% and 75% of the entire mass of the system. 
Although it is not our intention to discuss the dynamical evolution of this sys- 
tem it is worth mentioning that the core experiences a gravothermal runaway 
[46] at about t ~ 35, noticeable as the deepest point reached by the lower solid 
curve (5% Lagrangian radius), whereas the outer parts of the cluster continues 
to expand, mainly due to binary heating. 

A further analysis of simulation EqM using CWaT, is presented in fig. 7. For 
each star we compute the CWT for the duration of the simulations. In fig. 7 
we then show the fraction of stars (/eqm) which has a x > 0.1 (top curve), 0.2, 



10 



I ' ^ ' ^ ' ^ ' ^ ' 1 

20 40 60 80 100 



t 

Fig. 7. The evolution of /eqM for simulation model EqM, the fraction of stars with 
X > 0.1, 0.2, 0.4 and 0.6 (from top to bottom). 

0.4 and 0.6 (bottom curve). The various curves are rather smooth and flat. 
Some caution is well taken here, though, as near the edges (within the flrst 
and last ~ 20 time units) the CWaT becomes unreliable. The range between 
t = 20 and t = 80, where the method is accurate, we notice that the fraction 
of stars with any degree of chaoticity is roughly constant. It is interesting to 
note that ^ 20% of the stars have x > 0.6 and that 30-40% of the stars 
have ax < 0.2. In §3.2, by the way, we demonstrated that a value of x ~ 0.6 
is a strong indication that the underlying orbits change on a very short time 
scale. The core collapse around t ~ 35 is not clearly noticeable. Core collapse 
has apparently no direct influence in the number of irregular orbits in the 
equal-mass iV-body system. 

We run a second simulation with the same number of particles, distributed 
initially using the same 3 dimensional density profile, in which we assign ran- 
domly masses to each of the stars from a power-law distribution. Wc opted for 
a [47] mass function which has the form of N{m) oc m~^"^^. Here N{m) is the 
number of stars formed per unit mass range. We apply this mass function over 
a range of two orders of magnitude, i.e, the most massive star is a hundred 
times more massive than the least massive stars. We call this model S. 

We perform several simulations using this mass function but with new initial 

realizations of the system. The run-to-run variations can be quite large, mainly 
due to the small number of stars. A consistent trend, however, is evident. 

The dotted curves in flg. 6 present the evolution of the fundamental radii 
containing 5%, 50% and 75% of the mass of simulation model S. The evolution 

of the structure of model S is quite different than for the simulation with equal 
masses (model EqM). In the former case core collapse occurs much earlier (at 
about t = 14), which is consistent with simulations of [48]. However, model S 
also shows a second core collapse near t ~ 63. Noticeable is the much more 
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Fig. 8. The evolution of the ratio oi q = /s//EqM- The sohd and dotted curves 
are for two different reahzations of the initial conditions of the iV-body system 
with initial parameters according to model S. From top to bottom the curves are 
computed by adopting x > 0.1 (top solid and dotted curves), 0.2, 0.4 and x > 0.6 
(bottom sold curve). To prevent clutter the latter value for x is presented only for 
one of the simulations. 



dramatic expansion of the outer parts of the cluster, which is driven by the 
heating of the cluster core by close bound pairs which form during the core 
collapse (exempli gratia see [42,49]). 

We compute the CWT for each star in model S, for which we show the results 
of two simulations. For a better comparison we present the ratio gs/EqM = 
/s//EqM in fig- 8. Here /s and /eqm are computed for the same value of x, 
being either 0.1, 0.2, 0.4 or 0.6. Here the bottom curve represents the ratio in 
the number of stars for x > 0.6, the two curves give gs/EqM for x > 0.1. The 
ratio qs/EqM < 1 for all values of x > 0.1 and over the entire time span of the 
simulations. 

The result of two simulations with a mass spectrum (solid and dotted curves, 
in fig. 8 show the same behaviour. The structure of the curves for x > 0.4 and 
X > 0.6 are also quite interesting. They continue to drop at a constant rate 
until model S reaches a second core collapse around t — 60-70. 

The comparison carried out here, though not the main objective of this pa- 
per, is quite interesting from an astrophysical perspective. As the multi-mass 
simulation systematically appears to have fewer stars on irregular orbits. This 
deficiency can be as small as g = 0.05 during core collapse (or there about). 
We are currently extending this analysis to larger iV-body systems on order 
to study core collapse in large A^-body simulations. 
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4 Conclusions 



We present a qualitative diagnostic to analyze irregular behaviour in particle 
based gravitational A^-body simulations. The method, dubbed CWaT, is based 
on continuous wavelet transforms and provides a direct and instantaneous 
diagnostic for irregular behaviour. 

In our analysis we demonstrate that CWaT is suitable for analyzing self grav- 
itating A'"-body simulations, as it does not require a long time series for the 
analysis and it provides an instantaneous indicator for the degree of chaos, 
allowing a quantification of the number of irregular orbits at any moment in 
time in a large A'"-body simulation. 

Wc applied CWaT successfully to analyse gravitational 3-body interactions. 
The results of the CWaT method for 3-body systems compares well with 
Poincare sections. Here CWaT provides a quahfication of the irregularity of 
the orbit, whereas for Poincare sections such an analysis is technically much 
harder to perform. 

Eventually we applied the CWaT method to several N-body systems with 100 
particles in which all stars had the same mass, and another set of simulations 
in which a mass spectrum was adopted. The multi mass system systematically 
produced fewer stars on irregular orbits, compared to the equal mass system. 
In addition, we noticed that during a deep core collapse the number of irregular 
orbits dropped even further, whereas in the relatively low density inter-core 
collapse stages the number of irregular orbits tends to increases again. 
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Appendix: Initial conditions 

The calculations for testing the method are performed with the kira integra- 
tor which is part of the starlab gravitational A^-body environment [37] (see 
also http : / /www .manybody . org/manybody/starlab .html"). The initial con- 
ditions used for the 3-body encounters used as examples in sections 3.1 and 
3.2 are presented in the tables 1 and 2, respectively. The numbers in the ta- 
bles are in standard A^-body units [43] . Exact initial conditions for each of the 
three encounters studied is available via http://modesta.science.uva.nl 
The CWaT source code is also upon request. 

Table .1 

Initial conditions of the non-resonant preserving 3-body interaction. All numbers 
are in A^-body units The initial orbital parameters are: mi = m2 = ra^ = 0.5, a = 
1, e = 0, 6 ~ 1.34, Voo ~ 1.24t'c where the critical velocity Vc corresponds to a kinetic 
energy equal to the binding energy of the binary. [43] . 



star 


mass 


x-pos 


y-pos 


z-pos 


x-vel 


y-vel 




1 


0.5 


-32.1139669 


-19.6019186 


2.424416360 


-0.170641 


0.431148 


-0.0226399 


2 


0.5 


-32.6179535 


-20.465630 


2.424416360 


0.665287 


-0.117690 


-0.0226399 


3 


0.5 


64.73192040 


40.067549 


-4.84883272 


-0.494646 


-0.313458 


0.0452799 



Table .2 

Initial conditions of the democratic resonant preservation 3-body interaction. The 
initial encounter parameters are mi = m2 = = 0.5, a = 1, e = 0, 6 ~ 0.28, foo — 
0.26uc- 



star 


mass 


x-pos 


y-pos 


z-pos 


x-vel 


y-vel 


z-vel 


1 


0.5 


11.02969588 


-25.8888731 


25.3364950 


-0.498927 


0.206499 


-0.0524768 


2 


0.5 


10.74875682 


-26.8485987 


25.3364950 


0.453819 


-0.0972697 


-0.0524768 


3 


0.5 


-21.7784527 


52.7374718 


-50.6729901 


0.045107 


-0.10923 


0.104954 
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